/* 
 * File:   tersoffprob.hpp
 * Author: mikdebian
 *
 * Created on July 17, 2013, 2:03 AM
 */

#ifndef TERSOFFPROB_HPP
#define	TERSOFFPROB_HPP

#include <util/common/util.hpp>
#include <problems/optlib/objective.hpp>
#include <problems/boxcon/common/boxconproblem.hpp>
#include <tersoff/tersoffutils.hpp>
#include <lattice/lattice.hpp>
#include <tersoff/tersoffparams.hpp>

#define MAX_TERSOFF_PARAMS 10

/**
 * Objective reference values
 */
struct TersoffRefVal {
    /**
     * Total cohesive energy
     */
    double mE;
    /**
     * Bulk module
     */
    double mB;
    /**
     * Shear module
     */
    double mC;
    /**
     * Lattice constant
     */
    double mLC;
    /**
     * C44 module
     */
    double mC44;
    /**
     * Kleinman displacement
     */
    double mKlein;

};

class TersoffObjective : public Objective <double> {
public:

    /**
     * Options that control the objective
     */
    struct Options {
        /**
         * Include energy
         */
        static const unsigned int INC_ENERGY = 1;
        /**
         * Include Bulk
         */
        static const unsigned int INC_BULK = 1 << 1;
        /**
         * Include shear
         */
        static const unsigned int INC_SHEAR = 1 << 2;
        /**
         * Include lattice constant
         */
        static const unsigned int INC_LATTICE_CONST = 1 << 3;
        /**
         * Include C44
         */
        static const unsigned int INC_C44 = 1 << 4;
        /**
         * Include Kleinman internal displacement
         */
        static const unsigned int INC_KLEINMAN_DISP = 1 << 5;
    };

    /**
     * Constructor
     * @param tutil tersoff utility
     * @param rval reference values
     * @param n number of variables
     * 
     */
    TersoffObjective(TersoffUtils* tutil, TersoffRefVal* rval, int n) : mTut(tutil), mRval(rval) {
        mOptions = 0;
        mOptions |= Options::INC_ENERGY;
        mOptions |= Options::INC_BULK;
        mOptions |= Options::INC_SHEAR;
        mOptions |= Options::INC_LATTICE_CONST;
        mOptions |= Options::INC_C44;
        mOptions |= Options::INC_KLEINMAN_DISP;
        Objective<double>::setDim(n);
        eps = 0.00001;
    }

    double func(const double* x) {
        double rv = 0, v;
        int n = Objective<double>::mDim;
        auto dev = [](double v, double rv) { double u = (v - rv)/rv; return u * u;};
        copyTersoffParams(x);
//        printf("rv = %lf\n", rv);
        if (mOptions & Options::INC_ENERGY) {
            v = mTut->energyPerAtom();
  //          printf("v = %lf, ref = %lf\n", v, mRval->mE);
            rv += dev(v, mRval->mE);
        }
    //    printf("rv = %lf\n", rv);
        if (mOptions & Options::INC_BULK) {
            v = mTut->computeBulk();
            rv += dev(v, mRval->mB);
        }
      //  printf("rv = %lf\n", rv);
        if (mOptions & Options::INC_SHEAR) {
            v = mTut->computeShear();
            rv += dev(v, mRval->mC); ;
        }
        //printf("rv = %lf\n", rv);
        if (mOptions & Options::INC_LATTICE_CONST) {
            v = mTut->computeLatticeConst();
            rv += dev(v, mRval->mLC);
        }
       // printf("rv = %lf\n", rv);
        if (mOptions & Options::INC_C44) {
            double kli;
            v = mTut->computeC44Base(kli);
            rv += dev(v, mRval->mC44);
            if (mOptions & Options::INC_KLEINMAN_DISP) {
                rv += dev(kli, mRval->mKlein);
            }
        }
        //printf("v = %lf,rv = %lf\n", v, rv);
        return rv;
    }

    void grad(const double* x, double* g) {
        int n = Objective<double>::mDim;
        double y[n];
        VecUtils::vecCopy<double>(n, x, y);
        double fx0 = func(x);
        for (int i = 0; i < n; i++) {
            y[i] += eps;
            g[i] = (func(y) - fx0) / eps;
            y[i] -= eps;
        }
    }

    /**
     * Setup scale
     * @param s scale
     */
    void setupScale(double* s) {
        int n = Objective<double>::mDim;
        for (int i = 0; i < n; i++) {
            mScale[i] = BNBABS(s[i]);
        }
        printf("scale: ");
        VecUtils::vecPrint(n, mScale);
    }

    /**
     * Scale the params
     * @param x params to scale
     */
    void scale(double* x) {
        int n = Objective<double>::mDim;
        for (int i = 0; i < n; i++) {
            x[i] *= mScale[i];
        }
    }

    /**
     * Rescale the vector
     * @param x vector to rescale
     */
    void rescale(double* x) {
        int n = Objective<double>::mDim;
        for (int i = 0; i < n; i++) {
            x[i] /= mScale[i];
        }
    }

    /**
     * Returns options
     * @return options
     */
    unsigned int& getOptions() {
        return mOptions;
    }

private:

    void copyTersoffParams(const double* x) {
        int n = Objective<double>::mDim;
        double* t = (double*) mTut->getTpar();
        for (int i = 0; i < n; i++) {
            t[i] = mScale[i] * x[i];
            //  printf("x[%d] = %lf, mScale[%d] = %lf, t[%d] = %lf\n", i, x[i], i, mScale[i], i, t[i]);
        }
        //VecUtils::vecPrint(n, t);
    }

    unsigned int mOptions;
    TersoffUtils* mTut;
    TersoffRefVal* mRval;
    double eps;
    double mScale[MAX_TERSOFF_PARAMS];
};

struct TersoffProblem : public BoxconProblem <double> {
public:
    typedef BoxconProblem <double> BASE;

    TersoffProblem(TersoffObjective* obj, Box<double>& box) {
        BASE::mObj = obj;
        BASE::mBox = box;
    }

private:




};



#endif	/* TERSOFFOBJ_HPP */

